clear;clc;
disp('By Jun Luo')
s1=zeros(1,1);
a1=input('Please enter the lattice constant in the x direction for the first structure:');
b1=input('Please enter the lattice constant in the y direction for the first structure:');
xita=input('Please enter the angle between the x and y directions (0 < angle < 180):');
xit=xita*pi/180;
a=[cos(xit/2)*a1 sin(xit/2)*a1];
b=[cos(xit/2)*b1 sin(xit/2)*b1*(-1)];
N1=input('Please enter the number of atoms in the first structure:');
k=0;
v1=2;
for m1=-v1:v1
	for n1=-v1:v1
		for m2=-v1:v1
			for n2=-v1:v1
				k=k+1;
				A=m1.*a+n1.*b;
				B=m2.*a+n2.*b;
				s1(k,1)=m1;
				s1(k,2)=n1;
				s1(k,3)=m2;
				s1(k,4)=n2;
				s1(k,5)=norm(A);
				s1(k,6)=norm(B);
				s1(k,7)=(acos(dot(A,B)/(norm(A)*norm(B))))*180/pi;
				l1=real(s1);
			end
		end
	end
end

s2=zeros(1,1);
a2=input('Please enter the lattice constant in the x direction for the second structure:');
b2=input('Please enter the lattice constant in the y direction for the second structure:');
beita=input('Please enter the angle between the x and y directions (0 < angle < 180):');
beit=beita*pi/180;
c=[cos(beit/2)*a2 sin(beit/2)*a2];
d=[cos(beit/2)*b2 sin(beit/2)*b2*(-1)];
N2=input('Please enter the number of atoms in the second structure:');
k=0;
v2=2;
for m1=-v2:v2
	for n1=-v2:v2
		for m2=-v2:v2
			for n2=-v2:v2
				k=k+1;
				A=m1.*c+n1.*d;
				B=m2.*c+n2.*d;
				s2(k,1)=m1;
				s2(k,2)=n1;
				s2(k,3)=m2;
				s2(k,4)=n2;
				s2(k,5)=norm(A);
				s2(k,6)=norm(B);
				s2(k,7)=(acos(dot(A,B)/(norm(A)*norm(B))))*180/pi;
				l2=real(s2);
			end
		end
	end
end
mis_ang=input('Please enter the maximum angle mismatch you can accept (unit: deg):');
mis_rate=input('Please enter the maximum lattice mismatch rate you are allowed (generally 0.05):');


e1 = l1(:,7); e2 = l2(:,7);
t1 = l1(:,5); t2 = l2(:,5);
y1 = l1(:,6); y2 = l2(:,6);
s3=zeros(1,1);
k=0;
for m=1:(2*v1+1).^4
	for n=1:(2*v2+1).^4
		k=k+1;
		s3(k,1:4) = l1(m,1:4);
		s3(k,5:8) = l2(n,1:4);
		s3(k,9) = e1(m,1) - e2(n,1);
		s3(k,10) = e1(m,1);
		s3(k,11) = e2(n,1);
		sss1 = t1(m,1) - t2(n,1);
		sss2 = y1(m,1) - y2(n,1);
		s3(k,12) = sss1./t1(m,1);
		s3(k,13) = sss1./t2(n,1);
		s3(k,14) = sss2./y1(m,1);
		s3(k,15) = sss2./y2(n,1);
		k
	end
end

qw1=abs(s3(:,9));
qw2=abs(s3(:,12));
qw3=abs(s3(:,13));
qw4=abs(s3(:,14));
qw5=abs(s3(:,15));
index=find(qw1<mis_ang & qw2<mis_rate & qw3<mis_rate & qw4<mis_rate &qw5<mis_rate);
s4=s3(index,:);






k=0;
line_number1 = size(s4,1);
for o=1:line_number1    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%每次都要输入，这里1144是粗筛选过后s4的列数,
	if(s4(o,1)~=s4(o,3) | s4(o,2)~=s4(o,4) | s4(o,5)~=s4(o,7) | s4(o,6)~=s4(o,8));   %有一些夹角为0，因为配对了至少一对同样的晶格矢量（即相当于MS中输入了ab和ab），这里把0度的pass，但是180的没有pass
		k=k+1;
		s5(k,:)=s4(o,:);
	else
		llo=0;
	end
end

%如果这里算完，筛选出来的结果还是过多，必须要从原子的个数角度再进行筛选
aa1=[cos(xit/2)*a1 sin(xit/2)*a1 0];  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%每次都要输入，定义第一个材料，在没有变化晶格矢量时。注意：这里都是默认c方向垂直于ab面
bb1=[cos(xit/2)*b1 sin(xit/2)*b1*(-1) 0]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%每次都要输入
S1=norm(cross(aa1,bb1));  %叉乘法直接算第一个表面模型的面积（晶格矢量没有变化时）


aa2=[cos(beit/2)*a2 sin(beit/2)*a2 0];  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%每次都要输入，定义第二个材料，在没有变化晶格矢量时。注意：这里都是默认c方向垂直于ab面
bb2=[cos(beit/2)*b2 sin(beit/2)*b2*(-1) 0];  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%每次都要输入
S2=norm(cross(aa2,bb2));  %叉乘法直接算第二个表面模型的面积（晶格矢量没有变化时）

line_number2 = size(s5,1);
for i=1:line_number2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%每次都要输入，988为s5的列数
	m1 = s5(i,1); n1 = s5(i,2); m2 = s5(i,3); n2 = s5(i,4);
	AA1 = m1.*aa1+n1.*bb1;
	BB1 = m2.*aa1+n2.*bb1;
	S1_v1 = cross(AA1,BB1);
	S1_v = norm(S1_v1);
	N1_v = S1_v./S1.*N1;
	s5(i,16) = N1_v;    %材料1改变晶格矢量后原子的个数


	j1=s5(i,5); k1=s5(i,6); j2=s5(i,7); k2=s5(i,8);
	AA2 = j1.*aa2+k1.*bb2;
	BB2 = j2.*aa2+k2.*bb2;
	S2_v1 = cross(AA2,BB2);
	S2_v = norm(S2_v1);
	N2_v = S2_v./S2.*N2;
	s5(i,17)=N2_v;    %材料2改变晶格矢量后原子的个数
	s5(i,18)=N1_v+N2_v;
end
[~,index]=sort(s5(:,18));  %按照原子个数排序
s6=s5(index,:);

k=0;                    %最后的筛选，新晶格矢量的配对是纯数学的两两组合，会遇到组合后两两平行，此时原子数为0，没有物理意义，将其剔除
line_number3 = size(s6,1);
for o=1:line_number3
	if(s6(o,16)~=0 & s6(o,17)~=0);
		k=k+1;
		s7(k,:)=s6(o,:);
	else
		llkjo=0;
	end
end

save heterostructure.txt -ascii s7;
